GLMM example 9
Murray Logan
14 August, 2021
Preparations
Load the necessary libraries
library(car) #for regression diagnostics
library(broom) #for tidy output
library(broom.mixed) #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot) #for outputs
library(knitr) #for kable
library(effects) #for partial effects plots
library(ggeffects) #for effects plots in ggplotjk
library(emmeans) #for estimating marginal means
library(MASS) #for glm.nb
library(MuMIn) #for AICc
library(tidyverse) #for data wrangling
library(DHARMa) #for assessing dispersion etc
library(glmmTMB) #for glmmTMB
library(performance) #for diagnostic plots
library(see) #for diagnostic plots
library(ordinal) #for ordinal models
theme_set(theme_classic())Scenario
Read in the data
hughes <- read_csv('../data/hughes.csv', trim_ws=TRUE) %>%
janitor::clean_names()#bleaching data for 2016## Parsed with column specification:
## cols(
## REEF = col_character(),
## HABITAT = col_character(),
## SECTOR = col_character(),
## SCORE = col_double()
## )glimpse(hughes)## Rows: 3,394
## Columns: 4
## $ reef <chr> "09-357", "09-357", "09-357", "09-357", "09-357", "09-366", "…
## $ habitat <chr> "C", "C", "F", "F", "U", "C", "L", "C", "C", "L", "C", "C", "…
## $ sector <chr> "North", "North", "North", "North", "North", "North", "North"…
## $ score <dbl> 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 2…| reef | habitat | sector | score |
|---|---|---|---|
| 09-357 | C | North | 4 |
| 09-357 | F | North | 4 |
| 09-357 | U | North | 3 |
Data preparation
hughes <- hughes %>%
mutate(oscore = factor(score, ordered = TRUE),
habitat = factor(habitat),
sector = factor(sector, levels = c('North', 'Central', 'South')),
reef = factor(reef))For a categorical variable, it will apply a treatment contrast.
hughes.sum <- hughes %>%
count(sector, habitat, oscore) %>%
group_by(sector, habitat) %>%
mutate(prop = prop.table(n),
oscore = factor(oscore, levels = rev(levels(oscore))))Exploratory Data Analysis
ggplot(data = hughes.sum, aes(y = prop, x = habitat)) +
geom_bar(stat = 'Identity', aes(fill = oscore), color = 'black') + #'Identity' means don't do anything on it, multiply by 1
facet_grid(~sector) +
theme_bw() +
theme(panel.spacing.y = unit(10, 'pt'))
Fit models
hughes.clmm <- ordinal::clmm(oscore ~ habitat*sector + (1|reef),
data = hughes)
#partial plot
plot_model(hughes.clmm, type = 'eff',
terms = c('sector', 'habitat'))
#summary
summary(hughes.clmm)## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: oscore ~ habitat * sector + (1 | reef)
## data: hughes
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 3394 -3399.83 6831.67 1692(14086) 2.40e-03 6.5e+02
##
## Random effects:
## Groups Name Variance Std.Dev.
## reef (Intercept) 7.812 2.795
## Number of groups: reef 809
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## habitatF -0.3714 0.1930 -1.924 0.054304 .
## habitatL -1.4684 0.1909 -7.693 1.44e-14 ***
## habitatU 1.5678 0.4732 3.313 0.000924 ***
## sectorCentral -2.2063 0.3630 -6.078 1.21e-09 ***
## sectorSouth -7.1817 0.3017 -23.807 < 2e-16 ***
## habitatF:sectorCentral -0.6649 0.3389 -1.962 0.049745 *
## habitatL:sectorCentral 0.3329 0.3397 0.980 0.327179
## habitatU:sectorCentral -1.0853 0.6595 -1.646 0.099831 .
## habitatF:sectorSouth -0.3648 0.2494 -1.463 0.143585
## habitatL:sectorSouth -0.4088 0.2541 -1.609 0.107574
## habitatU:sectorSouth -0.7203 0.5323 -1.353 0.176024
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 -9.1435 0.2725 -33.553
## 1|2 -5.8773 0.2291 -25.659
## 2|3 -3.1539 0.1971 -15.998
## 3|4 -0.7224 0.1784 -4.049exp(-9.1435) #the odds for habitat C for having 0 bleaching compared to having a higher bleaching that that is 0.0001069125%## [1] 0.0001069125exp(1.5678) #odds of having higher bleaching score than habitat C in North sector## [1] 4.7960851/exp(-7.1817)## [1] 1315.142emmeans(hughes.clmm, ~oscore|habitat+sector, mode = 'prob') #gives a weighed average so the sum of each column probability is 1## habitat = C, sector = North:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 1.07e-04 2.91e-05 Inf 4.98e-05 1.64e-04
## 2 2.69e-03 6.12e-04 Inf 1.49e-03 3.89e-03
## 3 3.81e-02 7.18e-03 Inf 2.41e-02 5.22e-02
## 4 2.86e-01 3.28e-02 Inf 2.22e-01 3.50e-01
## 5 6.73e-01 3.93e-02 Inf 5.96e-01 7.50e-01
##
## habitat = F, sector = North:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 1.55e-04 4.80e-05 Inf 6.10e-05 2.49e-04
## 2 3.89e-03 1.05e-03 Inf 1.83e-03 5.95e-03
## 3 5.42e-02 1.25e-02 Inf 2.97e-02 7.88e-02
## 4 3.55e-01 4.46e-02 Inf 2.67e-01 4.42e-01
## 5 5.87e-01 5.66e-02 Inf 4.76e-01 6.98e-01
##
## habitat = L, sector = North:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 4.64e-04 1.39e-04 Inf 1.91e-04 7.36e-04
## 2 1.16e-02 3.01e-03 Inf 5.67e-03 1.74e-02
## 3 1.44e-01 2.90e-02 Inf 8.75e-02 2.01e-01
## 4 5.22e-01 2.59e-02 Inf 4.71e-01 5.73e-01
## 5 3.22e-01 5.16e-02 Inf 2.21e-01 4.23e-01
##
## habitat = U, sector = North:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 2.23e-05 1.20e-05 Inf -1.20e-06 4.58e-05
## 2 5.62e-04 2.89e-04 Inf -4.40e-06 1.13e-03
## 3 8.24e-03 4.07e-03 Inf 2.70e-04 1.62e-02
## 4 8.31e-02 3.66e-02 Inf 1.13e-02 1.55e-01
## 5 9.08e-01 4.09e-02 Inf 8.28e-01 9.88e-01
##
## habitat = C, sector = Central:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 9.70e-04 3.47e-04 Inf 2.91e-04 1.65e-03
## 2 2.38e-02 7.70e-03 Inf 8.75e-03 3.89e-02
## 3 2.55e-01 5.64e-02 Inf 1.44e-01 3.65e-01
## 4 5.36e-01 2.33e-02 Inf 4.90e-01 5.82e-01
## 5 1.85e-01 4.81e-02 Inf 9.05e-02 2.79e-01
##
## habitat = F, sector = Central:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 2.73e-03 1.07e-03 Inf 6.35e-04 4.82e-03
## 2 6.42e-02 2.22e-02 Inf 2.08e-02 1.08e-01
## 3 4.55e-01 6.87e-02 Inf 3.21e-01 5.90e-01
## 4 4.03e-01 6.64e-02 Inf 2.73e-01 5.34e-01
## 5 7.45e-02 2.54e-02 Inf 2.47e-02 1.24e-01
##
## habitat = L, sector = Central:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 3.01e-03 1.24e-03 Inf 5.81e-04 5.45e-03
## 2 7.04e-02 2.55e-02 Inf 2.04e-02 1.20e-01
## 3 4.73e-01 7.06e-02 Inf 3.35e-01 6.12e-01
## 4 3.85e-01 7.20e-02 Inf 2.44e-01 5.26e-01
## 5 6.79e-02 2.49e-02 Inf 1.92e-02 1.17e-01
##
## habitat = U, sector = Central:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 5.99e-04 3.32e-04 Inf -5.08e-05 1.25e-03
## 2 1.49e-02 7.86e-03 Inf -5.35e-04 3.03e-02
## 3 1.78e-01 7.42e-02 Inf 3.21e-02 3.23e-01
## 4 5.38e-01 2.74e-02 Inf 4.84e-01 5.92e-01
## 5 2.69e-01 1.04e-01 Inf 6.57e-02 4.72e-01
##
## habitat = C, sector = South:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 1.23e-01 2.17e-02 Inf 8.08e-02 1.66e-01
## 2 6.63e-01 2.02e-02 Inf 6.24e-01 7.03e-01
## 3 1.96e-01 2.99e-02 Inf 1.37e-01 2.54e-01
## 4 1.59e-02 3.41e-03 Inf 9.26e-03 2.26e-02
## 5 1.56e-03 3.80e-04 Inf 8.18e-04 2.31e-03
##
## habitat = F, sector = South:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 2.27e-01 3.88e-02 Inf 1.51e-01 3.03e-01
## 2 6.58e-01 2.29e-02 Inf 6.13e-01 7.03e-01
## 3 1.07e-01 2.12e-02 Inf 6.50e-02 1.48e-01
## 4 7.71e-03 1.89e-03 Inf 4.00e-03 1.14e-02
## 5 7.49e-04 2.04e-04 Inf 3.50e-04 1.15e-03
##
## habitat = L, sector = South:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 4.79e-01 5.59e-02 Inf 3.69e-01 5.88e-01
## 2 4.81e-01 4.82e-02 Inf 3.87e-01 5.76e-01
## 3 3.71e-02 8.58e-03 Inf 2.03e-02 5.40e-02
## 4 2.48e-03 6.57e-04 Inf 1.19e-03 3.77e-03
## 5 2.40e-04 6.94e-05 Inf 1.03e-04 3.76e-04
##
## habitat = U, sector = South:
## oscore prob SE df asymp.LCL asymp.UCL
## 1 5.68e-02 1.59e-02 Inf 2.56e-02 8.81e-02
## 2 5.55e-01 5.32e-02 Inf 4.51e-01 6.60e-01
## 3 3.48e-01 5.68e-02 Inf 2.36e-01 4.59e-01
## 4 3.63e-02 1.03e-02 Inf 1.60e-02 5.65e-02
## 5 3.64e-03 1.15e-03 Inf 1.39e-03 5.89e-03
##
## Confidence level used: 0.95emmeans(hughes.clmm, ~habitat|sector, mode = 'mean.class') #mean bleaching for each sector, need to subtract 1 from everything because bleaching classes start from 0## sector = North:
## habitat mean.class SE df asymp.LCL asymp.UCL
## C 4.63 0.0467 Inf 4.54 4.72
## F 4.52 0.0702 Inf 4.39 4.66
## L 4.15 0.0846 Inf 3.99 4.32
## U 4.90 0.0455 Inf 4.81 4.99
##
## sector = Central:
## habitat mean.class SE df asymp.LCL asymp.UCL
## C 3.88 0.1186 Inf 3.65 4.11
## F 3.48 0.1385 Inf 3.21 3.75
## L 3.44 0.1470 Inf 3.16 3.73
## U 4.06 0.1933 Inf 3.68 4.44
##
## sector = South:
## habitat mean.class SE df asymp.LCL asymp.UCL
## C 2.11 0.0556 Inf 2.00 2.22
## F 1.90 0.0615 Inf 1.78 2.02
## L 1.56 0.0647 Inf 1.44 1.69
## U 2.37 0.0938 Inf 2.19 2.56
##
## Confidence level used: 0.95emmeans(hughes.clmm, ~habitat|sector, mode = 'mean.class') %>%
pairs() #to directly compare the habitats by sector## sector = North:
## contrast estimate SE df z.ratio p.value
## C - F 0.105 0.0571 Inf 1.838 0.2557
## C - L 0.476 0.0675 Inf 7.063 <.0001
## C - U -0.269 0.0555 Inf -4.849 <.0001
## F - L 0.372 0.0823 Inf 4.513 <.0001
## F - U -0.374 0.0748 Inf -5.005 <.0001
## L - U -0.746 0.0881 Inf -8.464 <.0001
##
## sector = Central:
## contrast estimate SE df z.ratio p.value
## C - F 0.397 0.1062 Inf 3.740 0.0011
## C - L 0.435 0.1076 Inf 4.045 0.0003
## C - U -0.180 0.1695 Inf -1.061 0.7134
## F - L 0.038 0.1360 Inf 0.279 0.9924
## F - U -0.577 0.1890 Inf -3.053 0.0122
## L - U -0.615 0.1916 Inf -3.210 0.0073
##
## sector = South:
## contrast estimate SE df z.ratio p.value
## C - F 0.212 0.0457 Inf 4.641 <.0001
## C - L 0.545 0.0488 Inf 11.173 <.0001
## C - U -0.265 0.0801 Inf -3.313 0.0051
## F - L 0.333 0.0557 Inf 5.981 <.0001
## F - U -0.477 0.0860 Inf -5.547 <.0001
## L - U -0.810 0.0861 Inf -9.412 <.0001
##
## P value adjustment: tukey method for comparing a family of 4 estimatesnewdata <- emmeans(hughes.clmm, ~habitat|sector,
mode = 'mean.class') %>%
as.data.frame() %>%
mutate(across(c(mean.class, asymp.LCL, asymp.UCL), function(x) x-1))
ggplot(newdata) +
geom_hline(yintercept = 1, linetype = 'dashed', size = 0.1) +
geom_hline(yintercept = 2, linetype = 'dashed', size = 0.1) +
geom_hline(yintercept = 3, linetype = 'dashed', size = 0.1) +
geom_pointrange(aes(y = mean.class, x = habitat,
ymin = asymp.LCL, ymax = asymp.UCL)) +
facet_grid(~sector) +
scale_y_continuous('Bleaching score', breaks = (0:4), labels = 0:4,
limits = c(0, 4), expand = c(0, 0)) +
theme_bw()